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Coalescing binary black holes are a primary science target of ground-based gravitational-wave 
detectors. Detecting gravitational waves from these binaries and understanding the properties of 
the waves' sources require detailed knowledge of the expected waveforms. This paper presents 
a catalog of numerical binary black-hole simulations that represents a major advance toward the 
application of numerical relativity to gravitational- wave data analysis. The 171 simulations in 
the catalog include 90 processing binaries, mass ratios up to 8:1, orbital eccentricities from a few 
percent to 10 -5 , and black-hole spins up to 97% of the theoretical maximum. The simulations 
follow more orbits (up to 34) than previous simulations, allowing a more reliable connection to 
approximate analytical waveforms, which are used to extend numerical waveforms to span the entire 
frequency range of a detector. The expanded parameter-space coverage of this catalog, together 
with the length and accuracy of the included waveforms, will facilitate the development of new (and 
improvement of existing) gravitational-wave template families, the study of detection efficiency of 
gravitational-wave searches, and the understanding of systematic biases in parameter estimation of 
detected gravitational waves. Formidable challenges remain; for example, precession complicates the 
connection of numerical and approximate analytical waveforms, and vast regions of the parameter 
space remain unexplored. 

PACS numbers: 04.25.D-, 04.25.dg, 04.30.-w, 04.30.Db 



INTRODUCTION 

Gravitational waves (GW) from coalescing compact- 
object binaries — neutron stars (NS) and stellar-mass 
black holes (BH) — are primary targets for the next gen- 
eration of GW detectors, such as Advanced LIGO, Virgo 
and KAGRA PHI]- Detecting GWs from compact-object 
binaries requires high-quality, accurate theoretical wave- 
form models for GW template banks. Similarly, mea- 
suring source properties of detected signals ( "parameter 
estimation") relies on theoretical waveform models used 
in Markov Chain Monte Carlo algorithms [5 . 

For widely separated binaries, post-Newtonian (PN) 
calculations [5] provide gravitational waveforms to good 
accuracy. However, numerical simulations of the full Ein- 
stein equations are needed for the late inspiral, merger, 
and ringdown portions of the binary evolution. Such 
simulations are particularly important for stellar-mass 
BH-BH and BH-NS systems: Late inspiral and merger 
occur near LIGO's most sensitive frequency range, and 
although BHs might have high spins [5] , some of the 
spin contributions to the PN waveforms are known only 
to lower expansion order than the non-spinning terms [9]— 



This paper focuses on binary black holes (BBH). Nu- 
merical simulations of BBHs became possible eight years 
ago [33], with tremendous progress since (e.g., [3il 155] ). 
For best utility to GW astronomy, such simulations must 
satisfy several conditions: (i) sufficient accuracy; (ii) as- 
trophysical relevance, particularly low orbital eccentric- 
ity [SniinZ]; (hi) sufficient length (i.e., number of orbits) 
to connect reliably to PN waveforms; (iv) sufficiently 
dense coverage of relevant regions of parameter space. 

The precise requirements for these conditions are appli- 
cation dependent. For instance, 10 orbits is sufficient for 
GW detection at BBH parameters that are straightfor- 
ward to simulate (mass ratio q < 4, dimcnsionless spins 
\ = S/M 2 < 0.7 aligned with the orbital angular mo- 
mentum) [38] . On the other hand, parameter estimation 
can benefit from well over 100 orbits [551441] . Require- 
ments become more stringent with more extreme BH spin 
and mass ratio [40 . For precessing binaries, analogous 
studies have not even been performed because of a lack 
of inspiral-merger-ringdown waveform models. 

Satisfying all conditions (i) to (iv) is so difficult that, 
despite advances in numerical methods, simulations have 
barely reached the minimal desired quality. For instance, 
the world-wide NINJA-2 collaboration [42] resulted in 40 
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FIG. 1. Number of simulations vs. mass ratio and number 
of orbits. We separately list non-spinning, aligned-spin and 
precessing BBH systems. 



waveforms with an average length of ~9 orbits, essen- 
tially covering only two one-dimensional subspaces of the 
aligned spin region of the parameter space. The NR-AR 
collaboration [43, 44 is collecting and evaluating the ac- 
curacy of a smaller number of waveforms from systems 
with somewhat more general parameters than NINJA-2. 
A recent study [J5J reports about 80 simulations cover- 
ing approximately 7-12 orbits (^25 of them representing 
precessing binaries) and ^170 simulations lasting a few 
orbits each. 

This paper represents a major advance in waveform 
length and parameter-space coverage. We follow our ear- 
lier approach [4l)H4*9"] of computing black- hole simulations 
of the highest quality: longer than previous simulations, 
with higher accuracy and very low eccentricity. Our cat- 
alog contains 171 simulations, with almost all covering 
more than 12 orbits and the longest covering 34 orbits 
(Fig. [I]). 90 simulations represent precessing binaries. 

The remainder of this paper summarizes our compu- 
tational techniques, describes the catalog, and discusses 
applications and challenges for future work. 



TECHNIQUES 

The simulations described here are computed using the 
Spectral Einstein Code (SpEC) [50]. This multi-domain 
pseudospectral code includes an elliptic solver [5JJ and a 
hyperbolic evolution module [4"6"H4"5I 1521458"] to construct 
initial data and evolve it. 

Initial data are constructed [59] [60] to satisfy the Ein- 
stein constraint equations with two black holes in quasi- 
equilibrium 61 65j. Initial separation depends on the 
desired length of the inspiral, whereas the orbital fre- 
quency and the radial velocities of the black holes are 
determined iteratively [66, 67 to achieve eccentricities of 
~10~ 4 . 

SpEC employs singularity excision using a computa- 
tional domain extending from excision boundaries just 
inside each apparent horizon to some large radius of or- 



der several hundred M, where M is the total mass of the 
system. (Here G = c=l, so that time, length, and mass 
have the same dimensions.) Our first-order representa- 
tion [53] of the generalized harmonic formulation [581470] 
of Einstein's equations enables us to choose pure-outflow 
excision boundaries, which need no boundary conditions. 

Our outer boundary conditions [53] [TTJ [72] reduce the 
influx of constraint violations 73-79J and undesired in- 
coming gravitational radiation |80[ 181] . Running sim- 
ulations through merger robustly requires damped har- 
monic gauge conditions [55J [52] [53] > grids conforming to 
the shapes of the horizons [HI [551 [53 EH] , and adaptive 
mesh refinement [56] . After a common apparent hori- 
zon forms, the evolution is continued on a new computa- 
tional domain with a single excision boundary lying just 
inside the common apparent horizon and confornring to 
the common horizon's shape [HI [57] . 

Gravitational radiation is extracted from the simula- 
tions by two independent methods: the Newman-Penrose 
curvature scalar >F 4 as described in [HI IM] and the 
Regge- Wheeler- Zerilli (RWZ) gravitational wave strain 
him [5IH56"] . These waveforms, extracted at a series of 
finite-radius coordinate spheres, are extrapolated to infi- 
nite distance from the source using techniques developed 
in [87] with modifications suggested in [88] to account for 
precession. 



CATALOG 

The catalog is summarized in Fig. [2] Simulations 
0001-0065 correspond to initial-data set S from [57] . 
Most of these simulations have BH spins of either or 
0.5, oriented in different directions. Simulations 0066- 
0114 evolve several mass and spin configurations at mul- 
tiple eccentricities between 10 -4 and ~0.06. Simulations 
0115-0146 correspond to initial-data set S$ from |67j : 
they have random mass ratios in the range 1 < q < 2, 
random spin directions, and random spin magnitudes in 
the range < \ < 0.5. Simulation 0147 is a super- 
kick 89 - 91j configuration, with anti-parallel black hole 
spins tangential to the orbital plane. 

Simulations 0148-0165 use superposed Kerr-Schild ini- 
tial data [55J [52], which enable BH spins larger than 
X ~ 0.93, a bound that arises from other initial data 
constructions [65, 93J. Simulations 0149-0160 represent 
a sequence of equal-mass, equal-spin configurations with 
spins parallel or anti-parallel to the orbital angular mo- 
mentum (presented in more detail in [94 ) 

Simulations 0166-0171 represent unequal-mass non- 
spinning simulations and equal-mass, equal aligned spin 
simulations already published in [46] [95] [96] . The first 
published inspiral-merger-ringdown SpEC simulation (an 
equal-mass, zero-spin BBH [15] lasting about 4000M) 
is not included in the catalog, because we did not out- 
put enough information to apply our present extrapola- 
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FIG. 2. Properties of all simulations in the catalog. From top to bottom: dimensionless initial spin magnitudes xa,b\ angles 
0a,b between the initial spin vectors and the initial orbital angular momentum; angles 4>a,b between the line segment connecting 
the centers of the black holes and the initial spin vectors projected into the initial orbital plane; mass ratio q = Ma/Mb', number 
of orbits before merger; initial eccentricity e (triangles indicate an upper bound on e). 



tion techniques; however, improved equal-mass, zero-spin 
simulations are included (0001 and 0002). 

Figure [3] plots the gravitational- wave polarizations h+ 
and h x emitted into a certain sky direction, chosen such 
that h x vanishes for non-precessing systems. 

Figure El highlights two precessing systems, (i) Simula- 
tion 0035 is a q — 3 precessing configuration following 31 
orbits. This simulation proceeds through about 1.5 pre- 
cession cycles: the normal to the instantaneous orbital 
plane traces out the orange precession cone with open- 
ing angle 35° , whereas the BH spin of the more massive 
black hole traces out a precession cone with an opening 
angle of 144°. (ii) Simulation 0165 is a q = 6 configura- 
tion with both black holes spinning in generic directions, 
and spin magnitudes \A — 0.9, \b = 0.3. The orbital 
plane changes by almost 90°, and the spin direction of 
the smaller black hole traces out a spiral motion. 

Figure [5] shows the waveforms for simulation 0035 in 
two sky directions, highlighting the strong dependence of 
the waveforms on the orientation of the BBH relative to 
the line of sight to Earth. This figure also presents a con- 
vergence test, showing differences in the waveforms com- 
puted using different numerical resolutions. These differ- 
ences decrease rapidly with increasing resolution, illus- 
trating the exponential convergence provided by spectral 
numerical methods. Compared to our previous equal- 
mass non-spinning simulation [48, 49 and the unequal- 
mass non-spinning simulations |46] (all instrumental for 



various GW data-analysis applications [4"2l 155"! IST1 [55] ). 
0035 achieves a smaller cumulative phase error (0.07 rad) 
even though it is three times longer. Out of the 171 
simulations in the catalog, 57 (36) simulations were re- 
peated for at least two (three) resolutions. We estimate 
the accuracy of configurations with only one resolution 
by comparing with nearby configurations in the param- 
eter space; a more detailed discussion of the accuracy of 
waveforms in our catalog will be presented in [99] . 



DISCUSSION 

This paper presents the most comprehensive catalog of 
high-quality BBH simulations to date, enabling studies 
to help maximize the impact of GW detectors and to in- 
crease our understanding of GW sources and dynamical, 
strongly curved spacetime: 

Accuracy of PN precession equations: PN theory pre- 
dicts how spin and orbital angular momenta precess in 
generic binaries EH EH EH HMD ES E3 ETJ EH [3U 
I100H102J . The simulations here are long enough for de- 
tailed comparisons at different points in parameter space. 
Accurate modeling of precession dynamics is also ex- 
pected to be a key step toward constructing precessing 
waveform models [103J. 

PN accuracy studies, extended to a larger region of pa- 
rameter space and to an earlier stage of the inspiral: Pre- 
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FIG. 3. Waveform polarizations h+ (blue) and h x (orange) in a sky direction parallel to the initial orbital plane of each 
simulation. All plots have the same horizontal scale, with each tick representing a time interval of 2000Af (equal to 0.2 s for a 
2OM BBH). 



vious studies consider only aligned spins [47l 11041 1105] 
and at best either 15-20 orbits [351 EH] or a few orbits at 
very large separation |106j . With the longer waveforms 
in this catalog, these studies can extend to the earlier 
part of the inspiral, where PN theory is expected to be 
more accurate, and to include a larger region of parame- 
ter space. 

Independent validation of existing analytical waveform 
models: Many waveform models [86l I107H111] are cali- 
brated against numerical relativity simulations — but usu- 
ally only with a small number of simulations, which are 
comparatively short (typically < 10 orbits). The new 
simulations here enable tests of these models at different 
points in parameter space and covering more cycles. 

Detection sensitivity: Following the approach taken in 
the NINJA projects [121 [SSI EH], our waveforms can 
be injected into GW detector noise to study the effi- 
ciency of GW data-analysis pipelines. Injections of pro- 
cessing and/or eccentric waveforms from this catalog can 
quantify the impact of precession and eccentricity on the 
detection sensitivity of current searches using circular, 
aligned-spin templates. The new waveforms will also help 
assess the performance of searches with precessing wave- 
form templates. 

Systematic errors in parameter estimation: Parame- 
ter estimation methods [5] currently use inspiral-only PN 
waveforms. Applying parameter estimation methods to 



the waveforms in this catalog will enable the systematic 
errors of this approach to be quantified. 

Construction of precessing inspiral-merger-ringdown 
waveform models: This catalog will enable the construc- 
tion of sophisticated phcnomenological or effective-one- 
body waveform models for precessing binaries, extending 
the first such effort |113j . 

Periastron advance can be studied in aligned spin bi- 
naries and generic binaries [1141 1115] . 

While this catalog will enable pioneering studies, ma- 
jor challenges remain for future work. First, for a wave- 
form to be most useful for data analysis, it must be con- 
nected to a PN waveform from the early inspiral, forming 
a hybrid waveform |42j that spans the entire frequency 
range of a GW detector. This is difficult for precessing 
configurations, because of both the complexity of pre- 
cessing PN waveforms and ambiguities in connecting PN 
binary parameters with the numerical binary parame- 
ters [55] . 

Second, most of the parameter space remains unex- 
plored. Only 24 configurations have mass ratio q > 3 
(cf . Fig. [I]) ; of these, only 5 are precessing, and almost 
none have a spinning smaller black hole. Spinning BBH 
systems for 5 < q < 10 are particularly interesting be- 
cause they serve as accurate proxies for BH-NS binaries. 
Furthermore, the catalog contains only three simulations 
(the only three to date [HI [M]) with spins \ > 0.93 (i.e., 




FIG. 4. Precessing simulations 0035 and 0165. Top: Pre- 
cession cones traced by the orbital angular momentum and 
the spin vectors of the BHs . The circles denote the initial 
location of the vectors. For 0035 (left), \b — 0. Bottom: 
Trajectories of the individual black holes, projected into the 
initial orbital plane. 



above the "Bowen-York limit" |116F 120j). 

Finally, some of the simulations here utilize several al- 
gorithmic improvements, most notably adaptive mesh re- 
finement (AMR) to dynamically adjust the number of 
pseudospectral collocation points. AMR achieves high 
accuracy but complicates convergence tests; therefore, 
some simulations do not display convergence properties 
as clearly as as in Fig. [5] This will be discussed within 
a more comprehensive analysis of this catalog in a forth- 
coming publication [99 . Instructions for how to access, 
understand, and interpret the catalog will be provided 
in [99]. 
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FIG. 5. Waveforms and a convergence test for configuration 
0035. Top, middle: waveforms h+ (blue) and h x (orange) 
for two different angles (8, <j>) of the detector relative to the 
initial orbital angular momentum. The left panels show the 
inspiral and the right panels show the merger on an expanded 
horizontal scale. Bottom: cumulative phase differences be- 
tween the dominant £ — 2, m = 2 spherical-harmonic mode 
of the waveform of a given numerical resolution (labeled LI 
through L5 in order of increasing resolution) compared to that 
of the highest resolution L6. The plot excludes the early part 
of the waveform where initial transients are present. 
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